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Abstract 

A discrete adjoint method is developed and demonstrated for aerodynamic design optimization on unstruc- 
tured grids. The governing equations are the three-dimensional Reynolds-averaged Navier-Stokes equations cou- 
pled with a one-equation turbulence model. A discussion of the numerical implementation of the flow and adjoint 
equations is presented. Both compressible and incompressible solvers are differentiated and the accuracy of the 
sensitivity derivatives is verified by comparing with gradients obtained using finite differences. Several simplify- 
ing approximations to the complete linearization of the residual are also presented, and the resulting accuracy of 
the derivatives is examined. Demonstration optimizations for both compressible and incompressible flows are 
given. 


Introduction 

As computational power has continued to advance in recent 
years, researchers have been able to extend the use of computa- 
tional tools to increasingly more complex problems. Computa- 
tional fluid dynamics (CFD) has been exploited as an analysis 
tool for some time, and is currently receiving attention as a de- 
sign optimization tool. Early attempts in CFD-based design prob- 
lems made use of finite-difference calculations to obtain sensitiv- 
ity information. This technique can be used to obtain the 
derivatives of all the flow quantities with respect to each design 
variable and can be easily retrofitted to existing flow solvers. 
One problem with this approach is the computational time re- 
quired. To obtain the design sensitivities for a system involving 
n design parameters using a central-difference approach re- 
quires well-converged solutions of In flow analysis problems. 
For complex cases with many design variables, this requirement 
may become prohibitive. Another problem often encountered 
with the finite-difference approach is the sensitivity of the deriv- 
atives to the choice of the step size. It is desirable to have a small 
step size so that the truncation error is minimal, while at the same 
time, avoid exceedingly small step sizes which could yield large 
cancellation errors. 

To mitigate the difficulties associated with the choice of step 
size used in finite differences, direct differentiation can be em- 
ployed. 14, 19, 25 In this approach, the sensitivity derivatives of all 
the variables in the flow field are obtained but the solution of a 
large linear system of equations for each design variable is re- 
quired. Therefore, for problems involving many design variables, 
obtaining the sensitivity derivatives can be expensive. 

In recent years, adjoint formulations have grown in popularity, 
and are rapidly being developed for use in aerodynamic design 
sensitivity computations (see e.g. Refs. 2, 4, 9, 10, 15, 18). The 
adjoint approach has the advantage of being able to compute cost 
function gradients at an expense independent of the number of 
design parameters. This feature makes adjoint methods ex- 
tremely attractive for problems involving a large number of de- 
sign variables. 
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In the current work, a discrete adjoint approach is used in an un- 
structured-grid framework to compute design sensitivities for prob- 
lems based on the Navier-Stokes equations. A one-equation turbu- 
lence model is used which is fully differentiated and coupled into 
the solution of the adjoint equations so that the resulting sensitivity 
derivatives are consistent with those obtained using finite differ- 
ences. In addition to a compressible solver, an incompressible for- 
mulation is also differentiated in order to accommodate a wide 
range of applications. The accuracy of the resulting derivatives is 
established and sample calculations are shown for two- and three- 
dimensional cases using both implementations. Conclusions and 
suggestions for future research are also given. 
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Mach number 

Outward-pointing normal to control volume 
Prandtl number 
Turbulent Prandtl number 
Pressure 

Vector of conserved variables 
Components of heat flux 
Residual for a control volume 
Reynolds number 
Function for turbulence model 
Magnitude of vorticity 
Temperature 
Time 

Cartesian components of velocity 
Volume of control volume 
Grid-point locations 
Cartesian coordinate directions 


Angle of attack 
Ratio of specific heats 
Boundary of control volume 
Constant for turbulence model 
Vector of Lagrange multipliers 
Laminar viscosity 
Turbulent viscosity 
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Dependent variable for turbulence model 
Density 

Constant for turbulence model 
Components of shear stress 
Function for turbulence model 
Weights for lift and drag in cost function 


Subscripts 

00 Free-stream quantity 

Governing Equations 

The governing equations are the three-dimensional Reynolds-av- 
eraged Navier-Stokes equations. In the present work, both the com- 
pressible and incompressible forms of these equations are consid- 
ered. For turbulent flows, the one-equation turbulence model of 
Spalart and Allmaras 27 is used. The compressible flow equations, 
as well as the equations for the turbulence model, are given in Ap- 
pendix A. 


Adjoint and Design Equations 

In the adjoint approach for design optimization, a cost function is 
defined and augmented with the flow equations as constraints to 
form a Lagrangian given by 

L(D, Q, X, A) = f(D,Q,X)+A T R(D,Q,X) (1) 

Here, f(D, Q, X) is the cost function to be minimized and I) is a 
vector of design variables. The vector of Lagrange multipliers (also 
known as costate variables) is denoted by A , and R is the residual 
of the discretized steady-state flow equations. The vector Q is the 


conserved variables and X represents the computational grid. Al- 
though not explicitly denoted in Eq. 1 , both Q and X are functions 
of the design variables. 

Differentiating Eq. 1 with respect to the design variables yields 

§={® + SA} + 

Because A is arbitrary, the terms multiplied by may be 

eliminated using the following equation L J 
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Equation 3 is a linear system which represents the discrete adjoint 
equation for the optimization problem. After the flow equations 
have been solved for Q , the adjoint equation can be solved for the 
unknown vector of Lagrange multipliers A . The remaining terms 
in Eq. 2 can be used to evaluate the sensitivity derivatives as fol- 
lows: 
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dD [az> IdD] axj [IdD] Laz>J Idxj J 

After the solution for the costate variables is obtained using Eq. 3, 
the vector containing all of the desired sensitivities can be evalu- 
ated as a single matrix-vector product, given by Eq. 4. 

Numerical Implementation 


Flow Equations 

The flow solvers used in the current work are described at length 
in Refs. 1,3, and 5. The codes use an implicit, upwind, finite- vol- 
ume discretization, in which the dependent variables are stored at 
the mesh vertices. Inviscid fluxes at cell interfaces are computed 
using the upwind schemes of Roe, 21 van Leer, 29 or Osher. 6 Viscous 
fluxes are formed using an approach equivalent to a central-differ- 
ence Galerkin procedure. Temporal discretization is performed 
using a backward- Euler time-stepping scheme, and multigrid accel- 
eration can be used for the two-dimensional codes. 5 

An approximate solution of the linear system of equations 
formed at each time step is obtained using several iterations of a 
point-iterative scheme in which the nodes are updated in an even- 
odd fashion, resulting in a Gauss- Seidel-type method. 

The turbulence model is solved separately from the flow equa- 
tions at each time step, using a backward-Euler time-stepping 
scheme. The resulting linear system is solved using the same point- 
iterative scheme employed for the flow equations. The turbulence 
model is integrated all the way to the wall without the use of wall 
functions. 

The incompressible solvers are based on an artificial compress- 
ibility formulation, and are described in Ref. 3. Time integration 
and solution of the linear system at each time step are performed in 
the same manner as described above. 

Adjoint and Design Equations 

The adjoint equation given in Eq. 3 represents a linear set of 
equations for the co state variables A . Although this system can be 
solved directly using GMRES, 22 a time-like derivative is added and 
the solution is obtained by marching in time, much like the flow 
solver: 
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where 


A" + 1 = A" + A" A (6) 

The time term can be used to increase the diagonal dominance for 
cases in which GMRES alone would tend to stall. This ultimately 
results in a more robust adjoint solver. 

Due to the large amount of code resulting from the linearization 
of the viscous terms and the turbulence model, these contributions 
are stored in the present implementation. Because the stencil for the 
inviscid contributions is larger, the linearization of these terms is 
recomputed at each step to avoid the need for extra storage and data 
structure. 

To precondition the linear system, an incomplete LU decomposi- 
tion of the matrix obtained from a first-order accurate discretization 
is used. The preconditioning is applied on the left and no fill-in is 
allowed (ILU[0]). 13 Numerical experiments using this precondi- 
tioner have shown that some turbulent cases are slow to converge. 
An alternate means of preconditioning that has often been found 
useful is to employ a point-iterative scheme similar to that used for 
the flow equations. This technique allows for continual improve- 
ment in preconditioning the first-order system but is most effective 
when the time step is small and the matrix is diagonally dominant. 

In the present work, the differentiation of both the flow equations 
and the turbulence model is accomplished by ’’hand differentiating" 
the code. For obtaining the solution of the adjoint equations, the 
turbulence model is tightly coupled during the solution process, 
whereas it is solved separately during the flow analysis. During de- 
velopment, various treatments of the turbulence model have been 
studied and it has been found that the close coupling of the turbu- 
lence model is required in order to obtain sensitivity derivatives 
consistent with those obtained using finite differences. This will be 
illustrated in a subsequent section. 

Cost Functions 

For both two and three dimensions, the cost function is com- 
posed of a linear combination of the lift and drag coefficients: 

/ = - c;f + co 2 (Q - C' d i (7) 


For two-dimensional calculations, a target pressure distribution can 
also be specified. 

The drag can be minimized while maintaining a specified lift by 
adjusting the weights associated with each term in Eq. 7. The 
weights must be chosen such that neither term dominates the other. 
The current method for choosing the initial weights is to simply set 
the ratio of co 2 to to be equal to the ratio of the lift to the drag: 
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During the design process, these weights may require adjustment. 
However, this avoids the need to solve separate adjoint equations 
for lift and drag. 

Design Variables 

For both incompressible and compressible flows, the angle of at- 
tack can be utilized as a design variable in addition to the shape. 
For compressible flows, the free-stream Mach number can also be 
specified as a design parameter. When the shape is evolving, the 
surfaces are parameterized as described in the following sections. 


Two-Dimensional Parameterization 

For two-dimensional cases, the geometry is described using 
B-splines and the coordinates of the control points are used as de- 
sign variables as described in Refs. 2 and 4. Translation and rota- 
tion of individual bodies can also be used as design variables. Al- 
though not discussed further here, a graphical interface has been 
developed which aids the user in selecting and placing limits on de- 
sign variables, as well as modifying target pressure distributions. 2 
This interface has proven to be very useful and helps reduce errors 
in setting up design cases. 

Th ree-Dimension al Parameterization 

The three-dimensional code has been coupled with a geometric 
parameterization scheme recently introduced by Samareh. 23 The 
method utilizes a free-form deformation technique similar to that 
used in the motion picture industry for animating digital images. 
Here, a Bezier net describing the changes in the geometry is placed 
around the baseline mesh. The control points in the net may be used 
directly as the design variables, or they may be further grouped into 
design variables such as camber, thickness, and twist. This parame- 
terization technique has been chosen for its ability to handle arbi- 
trary geometries and because the mesh generation process does not 
depend on a prior parameterization of the geometry. This allows 
meshes which have been previously generated solely for analysis to 
be utilized for design purposes. 

Geometric Constraints 

During the design process, the feasibility of the geometry is 
maintained by limiting the movement of the design variables. For 
the two-dimensional code, area and curvature constraints may also 
be placed on the geometry although these are not employed in the 
present work. The curvature constraints can be used to enforce a 
specified leading edge radius or to guarantee curvatures of a speci- 
fied sign. 

Grid Generation and Mesh Movement Strategy 

For all of the two-dimensional computations, the meshes have 
been generated using the method described in Ref. 17. In three di- 
mensions, the method of Ref. 20 is utilized. Both techniques em- 
ploy an advancing front methodology and generate good quality 
grids for both inviscid and viscous calculations. 

When the design process requires modifications to the surface 
geometry, the computational mesh must be deformed to reflect the 
changes. For inviscid flows, the mesh movement strategy is based 
on the spring analogy described in Ref. 30. The edges of the mesh 
are treated as tension- springs, and the following equation is solved 
using a Jacobi iteration process: 

^ K ij (Ax i - Axj) = 0 (9) 

je N t 

Here Ax t and Axj represent the change in the coordinates of nodes 
i and j from the initial mesh to the desired mesh. The spring con- 
stants K f j are assumed to be /^ 2 , where / is the length of the edge 
connecting node i to node j . Since this technique may result in 
crossed grid lines, the required shift of the surface coordinates is 
decomposed into a series of smaller movements (usually around 
10), and Eq. 9 is relaxed for each change in the surface. This strat- 
egy has been found to work well for Euler-based designs. 

For viscous meshes, the method described above is not adequate 
and can easily lead to crossing mesh lines and negative volumes. 
For these cases, the nodes near viscous surfaces are shifted by inter- 
polating the changes in the coordinates at the boundaries of the 
nearest surface triangle or edge. This technique is blended with a 
smoothing procedure so that away from the highly stretched cells 
near the surface, the mesh movement reverts to that of the proce- 
dure described above for inviscid meshes. Further details can be 
found in Ref. 4. 
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For two-dimensional viscous applications, the procedure de- 
scribed above is very robust and is capable of successfully deform- 
ing the mesh in response to large changes in the surface geometry. 
For multielement airfoils, there is a tendency to open "gaps" in the 
mesh between elements when the elements are allowed to translate 
away from one another. A similar problem occurs when the ele- 
ments move closer together in which case there is a "jamming" to- 
gether of mesh points. These difficulties are simply due to the fact 
that no additional mesh points are inserted or removed during the 
process so that as elements shift in relation to one another, voids 
can be created. This has not had a detrimental effect on the flow 
solver and can be remedied by periodically regenerating the mesh. 

For three-dimensional flows, the mesh movement procedure is 
inadequate when large changes in the geometry are required. In 
these cases, negative cell volumes can occur around the edges of 
the planform. Therefore, changes in the thickness and camber have 
been limited to only a few percent of the chord. Further research is 
required to develop a more reliable methodology for large geomet- 
ric changes in three dimensions. 

As the surface is deformed during the design process, there is a 
corresponding change in the interior mesh points as well. The effect 
of the changing grid is reflected through the mesh sensitivity terms 
given by dX/dD in Eq. 4. The computation of these terms is 
achieved by differentiating the mesh movement process described 
above. 

Optimization Technique 

The optimization technique used in all of the results presented 
below is the quasi -Newton method of Davidon- Fletcher- Powell. 8,11 
The current implementation of this technique, referred to as 
KSOPT, allows for multipoint optimization as well as both equality 
and inequality constraints. 32 For the present work, the multipoint 
capability is not utilized although this is an obvious future require- 
ment. 

Results and Discussion 

Consistency of Linearization 

During code development, great care has been taken at each step 
to ensure that the derivatives are consistent with those obtained 
using finite differences. In this section, the accuracy of the resulting 
derivatives is verified for compressible flow in both two and three 
dimensions. Similar results have been obtained for incompressible 
flows and are included in Appendix B. Comparisons are made be- 
tween derivatives computed using finite differences with those ob- 
tained using the adjoint method. When computing derivatives using 
finite differences, central-difference formulas are used with a step 
size of lxl 0~ , and the flow solver is converged to machine accu- 
racy. The two-dimensional results shown are calculated using the 
Osher flux function although similar accuracy is obtained using ei- 
ther flux-difference splitting or flux-vector splitting. For the three- 
dimensional linearizations, all results are obtained using the flux- 
difference splitting scheme of Roe. All of the results shown below 
are for turbulent flows although the consistency of derivatives has 
been verified for inviscid and laminar flows as well. 

Two-Dimensional Accuracy 

For demonstrating the consistency of the derivatives obtained 
using the adjoint formulation with those obtained using finite dif- 
ferences, two test cases are considered. The first case is a 2-element 
airfoil at a free-stream Mach number of 0.25, an angle of attack of 
1 ° , and a Reynolds number of 9 million based on the chord of the 
airfoil. The geometry has been chosen arbitrarily and is that given 
in Ref. 3 1 . 

The mesh used for this test has 4,901 nodes and is shown in Fig. 
1. The geometry of each airfoil is described with a third-order 
B-spline. The derivatives of the lift and drag coefficients with re- 
spect to the vertical and horizontal positions of four shape design 
variables have been obtained. The locations of the design variables 
are indicated by the solid circles shown in Fig. 2. As seen in the fig- 


ure, two of these design variables are located on the main airfoil 
and two are located on the flap. For each element, one design vari- 
able is located on the upper surface near the nose of the airfoil and 
one is located near the rear. A comparison of derivatives of the lift 
and drag coefficients with respect to changes in the vertical position 
of these design variables is shown in Tables 1 and 2, while Tables 3 
and 4 compare the derivatives of the lift and drag coefficient with 
respect to x- and y-translation of the flap. Note that this required 
two solutions of the adjoint equation — one for lift and one for 
drag. As seen, the derivatives obtained with the adjoint approach 
are in very good agreement with the finite-difference derivatives 
for all cases. Although not shown, similar accuracy is obtained for 
the derivatives with respect to horizontal changes in the control 
points. 



Figure 1. Mesh for 2-element airfoil used in assessment of two- 
dimensional design sensitivities. 
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Figure 2. Location of design variables for 2-element airfoil. 


Table 1. Accuracy of two-dimensional derivatives for lift coefficient 



Finite difference 

Adjoint 

% Diff. 

AC 

0.30965 

0.30962 

-0.010 

a 

0.12285 

0.12282 

-0.024 

Point A 

-1.0952 

-1.0952 

0.000 

Point B 

0.57480 

0.57480 

0.000 

Point C 

-2.1368 

-2.1366 

-0.009 

Point D 

0.76215 

0.76215 

0.000 


In order to further demonstrate the accuracy of the differentia- 
tion, a case containing transonic flow is examined. An RAE 2822 
airfoil is used at an angle of attack of 2.81° , a Mach number of 
0.75 , and a Reynolds number of 6.2 million. The mesh contains 
14,127 nodes and the spacing at the wall is 1x10” . The computed 
pressure distribution is shown in Fig. 3 along with the correspond- 
ing experimental data. 7 For this case, a strong shock is present on 
the upper surface which separates the flow immediately down- 
stream. The locations of the three design variables are shown by the 
filled circles in Fig. 4. The corresponding sensitivity derivatives for 
the lift coefficient with respect to a vertical movement of the con- 


4 


Table 2. Accuracy of two-dimensional derivatives for drag coefficient 



Finite difference 

Adjoint 

% Diff. 

AC 

-0.05029 

-0.05029 

0.000 

a 

0.00843 

0.00843 

0.000 

Point A 

0.21925 

0.21925 

0.000 

Point B 

-0.03489 

-0.03489 

0.000 

Point C 

0.17007 

0.17007 

0.000 

Point D 

0.06447 

0.06448 

0.016 


Table 3. Accuracy of two-dimensional derivatives of lift coefficient 
for flap translation 



Finite difference 

Adjoint 

% Diff. 

x-translation 

1.4226 

1.4232 

0.042 

y-translation 

-6.8991 

-6.8990 

-0.001 


Table 4. Accuracy of two-dimensional derivatives of drag 
coefficient for flap translation 



Finite difference 

Adjoint 

% Diff. 

x-translation 

0.02710 

0.02716 

0.221 

y-translation 

-0.24164 

-0.24163 

-0.004 



Figure 3. Pressure distribution for transonic RAE 2822 airfoil. 



Figure 4. Location of design variables for RAE 2822 airfoil. 


trol points are listed in Table 5 along with those for Mach number 
and angle of attack. The agreement with finite differences is very 
good, although the derivative for the design variable located in the 
separated region of the flow just downstream of the shock appears 
to be slightly inconsistent. However, numerical experiments using 
different step sizes have shown that finite differences for this con- 
trol point are somewhat sensitive to the perturbation level. For ex- 
ample, step sizes of 1 x 10~ 6 and 5x10 result in finite-differ- 
ence derivatives of 1.9187 and 1.9150 respectively. 


Table 5. Sensitivity derivatives for lilt coefficient for RAE 2822 airfoil 



Finite difference 

Adjoint 

% Diff. 

AC 

-3.0546 

-3.0546 

0.000 

a 

5.7614 

5.7615 

0.002 

Point A 

7.9826 

7.9814 

-0.015 

Point B 

1.9247 

1.9219 

-0.145 

Point C 

1.3283 

1.3283 

0.000 


Three-Dimensional Accuracy 

To verify the accuracy of the derivatives in three dimensions, a 
similar experiment is conducted. For this case, an ONERA M6 
wing 24 has been parameterized using 46 design variables describing 
the planform, twist, shear, thickness, and camber. The design vari- 
ables are depicted in Fig. 5 where twist and wing shear have been 
parameterized at five spanwise locations. The thickness and camber 
have also been parameterized using the six locations shown in the 
figure. The design variables describing the planform are not shown 
in the figure nor are thickness and camber design variables along 
the leading and trailing edges. The mesh used for these tests con- 
tains 16,391 nodes and 90,892 tetrahedra and is shown in Fig. 6. 
The flow conditions are an angle of attack of 2° , a Reynolds num- 
ber of 5 million based on the mean chord, and a Mach number of 
0.3 . In this test, the cost function is a combination of the lift and 
drag coefficients so that only one adjoint solution is required. The 
derivatives of the cost with respect to the angle of attack and the 
Mach number as well as the derivatives with respect to four of the 
shape parameterization variables are shown in Table 6. As can be 
seen, the consistency between the derivatives obtained with the ad- 
joint formulation and finite differences is excellent. Additional de- 
rivatives for the design variables depicted in Fig. 5 have also been 
verified with comparable accuracy. 


Twist 

Shear 



5 






















Table 7. Sensitivity derivatives for bit coefficient using various 
approximations (see Fig. 1) 



Figure 6. Grid used for assessment of three-dimensional design 
sensitivities. 


Table 6. Three-dimensional compressible derivatives 



Finite difference 

Adjoint 

% Diff. 

Mach Number 

0.00960 

0.00959 

-0.104 

a 

-0.03243 

-0.03243 

0.000 

Twist #3 

0.00965 

0.00965 

0.000 

Shear #3 

-0.04275 

-0.04277 

0.047 

Thickness #3 

-0.04011 

-0.04012 

0.025 

Camber #4 

-1.3174 

-1.3174 

0.000 


Linearization Approximations 

Due to the complexity in achieving accurate linearizations for 
use in Eqs. 3 and 4, one may consider the use of simplifying 
assumptions. Clearly, a great deal of effort can be avoided if certain 
terms may be neglected or replaced with simpler approximations 
without seriously compromising the accuracy of the results. The 
previous sections have established the accuracy of the derivatives 
obtained from the adjoint formulation using a consistent lineariza- 
tion of the flow solvers. This section will examine the accuracy of 
the derivatives obtained using several natural approximations. 
These numerical experiments are conducted in two dimensions 
using the test case and flow conditions used for Table 1. A discus- 
sion of each of the approximations is given below and some repre- 
sentative derivatives for vertical changes in the design variables are 
shown in Table 7. 

First-Order Adjoint Solution 

For second-order accurate schemes, the complete linearization of 
the inviscid contribution to the residual requires information from 
mesh points beyond the immediately adjacent nodes. This require- 
ment arises from having to form gradients of the dependent vari- 
ables at the nodes in order to extrapolate them to the faces of each 
control volume. This large stencil makes an exact linearization 
quite tedious. However, if the fluxes are formed using only nearest- 
neighbor information, the amount of coding required above the 
baseline flow solver is minimal. This corresponds to using a first- 
order accurate scheme for the convective terms, and may also result 
in a linear system that is easier to solve, as the bandwidth of the co- 
efficient matrix is reduced significantly. 



Finite 

difference 

Adjoint 

% Diff. 

Exact Linearization 




Point A 

-1.0952 

-1.0952 

0.000 

Point B 

0.57480 

0.57480 

0.000 

Point C 

-2.1368 

-2.1366 

-0.009 

Point D 

0.76215 

0.76215 

0.000 

x-translation (flap) 

1.4226 

1.4232 

0.042 

y- translation (flap) 

-6.8991 

-6.8990 

-0.001 

First-Order Adjoint 




Point A 

-1.0952 

0.34633 

-132 

Point B 

0.57480 

0.47104 

-18.1 

Point C 

-2.1368 

-0.10590 

-95.0 

Point D 

0.76215 

0.70373 

-7.67 

x-translation (flap) 

1.4226 

-0.68201 

-148 

y-translation (flap) 

-6.8991 

0.14801 

-102 

Constant p, 




Point A 

-1.0952 

-1.6844 

53.8 

Point B 

0.57480 

0.37262 

-35.2 

Point C 

-2.1368 

-2.2888 

7.11 

Point D 

0.76215 

0.57459 

-24.6 

x-translation (flap) 

1.4226 

1.1909 

-16.3 

y-translation (flap) 

-6.8991 

-7.6163 

10.4 


In Table 7, derivatives obtained using a first-order linearization 
of the convective terms are compared with those obtained from the 
linearization of the higher order residual. For these results, the first- 
order approximation is made in evaluating both Eqns. 3 and 4. 
Using this approximation, the derivative of the lift with respect to a 
vertical shift of the design variable towards the rear of the flap is 
within 8% of the correct value. However, the derivatives obtained 
by ignoring the higher order terms are generally highly inaccurate 
and several are of incorrect sign. The derivatives of incorrect sign 
would most certainly have an adverse effect on an optimization 
process, especially near a minimum. 

" Frozen " Turbulence Model 

An accurate linearization of the turbulence model can be diffi- 
cult to obtain. As seen from the equations given in Appendix A, 
there are many terms and additional functions that must be properly 
differentiated. These terms exhibit complex dependency on both 
the flow variables as well as the distance to the wall. By assuming 
that the turbulence model is "frozen", a significant reduction in the 
required level of effort may be obtained. This approach has been 
previously used in Refs. 15 and 26 for structured grid applications 
to airfoils and wings. In these references, successful optimizations 
have been performed although the accuracy of the derivatives has 
not been explicitly demonstrated. 

Results obtained by making the assumption of a constant eddy 
viscosity are listed in Table 7. While the computed sensitivities 
show a large amount of error when compared to finite differences, 
they all exhibit the correct sign. However, for derivatives 
associated with horizontal changes in these same design variables, 
several are of incorrect sign. For example, the finite-difference 
derivative obtained by perturbing point A in the horizontal 
direction is -0.18060 whereas the current approximations to the 
linearizations yield 0.30328. For this same design variable, the 
complete linearization yields -0.18052 which is less than 0.05 
percent in error. 

A similar technique that can be used to simplify the implementa- 
tion is to neglect the contributions from the turbulence model in 
Eq. 4. This is primarily motivated by the observation that the 
co state variable associated with the turbulence model is very small 
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and decays rapidly away from the body. Although not shown, nu- 
merical experiments indicate that the resulting accuracy is poor 
with many derivatives of incorrect sign. 

Extent of Mesh Sensitivities 

For each design variable, the evaluation of Eq. 4 requires a ma- 
trix-vector product of the co state variables with the linearization of 
the residual. This also includes computation of the mesh sensitivi- 
ties for each design variable. For large numbers of design variables 
and mesh points, this can potentially represent a significant expense 
due to the complicated linearization of the residual. Because nodes 
further away from the body are subjected to more moderate 
changes than those in the immediate vicinity of the surface, it may 
be possible to neglect terms in Eq. 4 that are sufficiently far from 
the body. This could help to reduce the cost of evaluating Eq. 4 by 
avoiding the need to include terms from every mesh point in the 
field. 

To investigate the validity of this assumption, a region around 
the surface of the airfoil is defined by first "tagging" the nodes on 
the surface and then identifying nodes that lie within a set number 
of grid layers adjacent to the surface. Outside of this region, the 
mesh sensitivities are set to zero to emulate the effect of neglecting 
all the contributions outside of the tagged region. The sensitivity 
derivatives for the lift coefficient with respect to vertical and hori- 
zontal translations of the flap are computed for a varying number of 
grid layers and the results are shown in Fig. 7. Here, Ci trans is the 
ratio of the approximate derivative to the derivative obtained by in- 
cluding the mesh sensitivities at every grid point in the domain. In 
this figure, the curve labeled n/n tot is the ratio of the number of 
nodes where mesh sensitivities are employed to the total number of 
nodes in the mesh. It should be noted that examining a single deriv- 
ative may not be representative of the behavior of the rest of the de- 
rivatives and an accurate computation of this derivative does not 
guarantee accuracy for the remaining derivatives. However, inaccu- 
racy of this derivative demonstrates that neglecting the full effects 
of the mesh sensitivities may have an adverse effect on other deriv- 
atives as well. 

As seen in the figure, the influence of the mesh sensitivities grad- 
ually decays away from the surface. Accurate results are obtained 
when the number of mesh layers is greater than approximately 15. 
At this point, about half the total number of points in the mesh is in- 
cluded in the layers so that a factor of 2 savings could be realized 
when evaluating Eq. 4. When many design variables are present, 
neglecting some of the mesh sensitivities could lead to a substantial 
savings in computer time. However, for the present study, the com- 
puter time required for evaluating Eq. 4 does not dominate the over- 
all optimization process so this strategy is not used. 



Design Examples 

Inviscid Multielement Airfoil 

The objective of the first example computation is to position the 
elements on a multielement airfoil to obtain a desired pressure dis- 
tribution. For this case, the target pressure distribution is obtained 
from analysis of a baseline configuration. 28 The individual ele- 
ments are then perturbed by translating in the x- and y-directions as 
well as by rotating by several degrees. The mesh used for this test 
contains 3,820 nodes with 1 93 nodes on the surface of the main ele- 
ment and 129 on each flap. The initial and target configurations are 
shown in Fig. 8. 

The initial pressure distributions over the elements are shown in 
Fig. 9 along with the target pressure distribution and the pressures 
obtained after 15 design cycles. A near-field view of the corre- 
sponding geometries are shown in Fig. 10. As seen in the figures, 
the target pressure distributions are obtained and each of the ele- 
ments returns to its desired position. This experiment has been suc- 
cessfully performed using both the compressible and incompress- 
ible solvers, although only results from the incompressible case are 
shown. 



Figure 8. Four-element airfoil in original and perturbed positions. 


Current 

Target 

Initial 




-2.1 -2.0 -1.9 -1.8 -1.7 -2-10 1 2 



x x 


Figure 9. Pressure distributions for 4-element airfoil. 


Figure 7. Extent of mesh sensitivity terms required for translation 
sensitivity accuracy. 
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Current 

Target 



-2.1 -2.0 -1.9 -1.8 -1.7 -2 




x 


x 


Figure 10. Geometries for 4-element airfoil. 


Turbulent Airfoil 

The objective for the next test case is to reduce the drag on the 
RAE 2822 airfoil. The initial flow conditions are a ffee-stream 
Mach number of 0.75 , an angle of attack of 2.81° , and a Reynolds 
number of 6.2 million based on the chord of the airfoil. These con- 
ditions correspond to those presented earlier for verifying the shape 
sensitivity derivatives. For this case, there are 47 active design vari- 
ables. The lift coefficient is held fixed at 0.7336 and the objective 
is to reduce the drag coefficient. After 10 design cycles, the drag 
has been reduced from 0.0263 to 0.0150 whereas after 20 design 
cycles, a modest improvement is further obtained, reducing the 
drag coefficient to 0.0149 . The initial and final pressure distribu- 
tions are shown in Fig. 1 1 and Mach number contours are shown in 
Fig. 12. As seen in the figures, the shock wave on the surface of the 
airfoil is eliminated although the curvature in the Mach contours for 
the final geometry indicate that a shock may form at off-design 
conditions. 


Multielement Airfoil 

For this case, the objective is to increase the downforce for a 
multielement airfoil used for open-wheel racing cars. This airfoil 
has been initially designed using the inviscid design techniques de- 
scribed in Ref. 12. The mesh has 15,446 nodes with 195 placed on 
the main element and 129 on the flap. The spacing at the wall is 
1x10” giving a y ~ 1 based on flat plate estimates. The Rey- 
nolds number is 2.4 million based on the chord of the airfoil and the 
angle of attack is held fixed at 12° which corresponds to the oper- 
ating point suggested in Ref. 12. After 25 cycles using 65 design 
variables, the downforce coefficient has increased from -2.3068 to 
-2.4379. The initial and final pressure distributions computed using 
the incompressible solver are shown in Fig. 13 along with the corre- 
sponding shapes. It is seen from the figure that the redesigned main 
element carries more downforce compared to the initial design, 
while the loading on the flap has decreased. Velocity contours and 
vectors, shown in Fig. 14, indicate that the region of separated flow 
on the rear of the flap has been reduced through the design process. 
Although not shown, comparable cases have been run using the 
compressible code with similar results. 



Figure 11. Initial and final pressure distributions for drag reduction 
on RAE 2822 airfoil. 
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Figure 12. Initial and final Mach contours for transonic airfoil 
optimization exercise. 




Figure 13. Initial and final pressure distributions for multielement 
airfoil. 




Figure 14. Velocity contours and vectors in flap region for 
multielement airfoil. 


Inviscid Drag Reduction for ONERA M6 Wing 

An example optimization is conducted for inviscid flow over the 
ONERA M 6 wing. 24 The tree-stream Mach number for this case is 
0.84 and the angle of attack is 3.06° . The mesh used for this com- 
putation consists of 53,961 nodes and 287,962 tetrahedrons. The 
contours for the initial and final density distribution on the surface 
of the wing are shown in Fig. 15 with the corresponding pressure 
distributions shown in Fig. 16. The objective of the optimization is 
to reduce the drag while maintaining a specified lift. For this de- 
sign, the angle of attack is allowed to change in addition to 10 
shape design variables (4 twist, 4 camber, and 2 thickness). The 
twist variables are located at the 4 outboard stations in Fig. 5 and 
allowed to increase or decrease by 1 ° . The thickness and camber 
variables at positions 3 and 4 are also design variables as is the 
camber at positions 5 and 6. Each of these is allowed to change by 2 
percent of the span. After 1 0 design cycles, the drag has been re- 
duced from 0.0182 to 0.0167 while the lift has been maintained. 
The pressure distribution shown in Fig. 16 indicates that the shock 
has weakened at all of the span wise stations. 




Figure 15. Initial and final density contours for inviscid wing 
design. 
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x/c x/c 



x/c x/c 


Figure 16. Initial and final pressure distributions for inviscid wing 
design. 

Turbulent ONERA M6 Wing Redesign 

A transonic wing design has been conducted using an ONERA 
M6 mesh consisting of 62,360 nodes and 355,814 tetrahedra. The 
mesh used for this study is extremely coarse and is not adequate for 
accurate computations; it serves merely as an initial demonstration 
for evaluating the methodology. The flow is assumed to be fully 
turbulent at a Mach number of 0.84, an angle of attack of 3.06° , 
and a Reynolds number of 5 million. For this flow field, a weak 
swept shock extends from the root leading edge and a normal shock 
is present further aft on the wing surface (see Fig. 17). The weak- 
ness of the shock is in large part due to the coarseness of the mesh; 
further refinement would yield a shock structure similar to that 
shown in the initial flow field in Fig. 15. The objective of the de- 
sign is to reduce the drag while holding the lift constant. For this 
case, thickness and camber have been allowed to vary at two chord- 
wise stations located at the mid-span of the wing. These design 
variables have been allowed to change up to 1 percent of the span 
of the wing. The angle of attack is also allowed to vary in order to 
maintain the original lift coefficient. 

After 10 design cycles, the drag coefficient is reduced from 
0.0200 to 0.0184. Density contours for the initial and final design 
are shown in Fig. 17. It is apparent from the increased spacing be- 
tween the contours that the strength of the shock at the midchord lo- 
cation is somewhat weaker for the final design which accounts for 
the lower drag. 

Summary and Concluding Remarks 

Compressible and incompressible versions of an unstructured 
mesh Navier-Stokes flow solver have been differentiated and the 
resulting derivatives have been verified by comparisons with finite 
differences in both two and three dimensions. In this implementa- 
tion, the turbulence model is fully coupled with the flow equations 
in order to achieve this consistency. The accuracy of a number of 
simplifying approximations to the linearizations of the residual 
have also been examined and none of the approximations yielded 
derivatives of acceptable accuracy and were often of incorrect sign. 

Efficient surface parameterizations have been utilized in both 
two and three dimensions, and the resulting codes have been inte- 
grated with an optimization package. Example optimizations have 
been demonstrated in both two and three dimensions. 




Figure 17. Initial and final density contours for wing. 

In order for large scale optimization to become routine, the bene- 
fits of parallel architectures should be exploited. Although the 
three-dimensional flow solver has been parallelized using compiler 
directives, the parallel efficiency is under 50 percent. Clearly, par- 
allel versions of the codes will have an immediate impact on the 
ability to design realistic configurations on fine meshes, and this ef- 
fort is currently underway (see e.g. Ref. 16). 

Another area that requires future work is the incorporation of 
multipoint optimization capability for designing geometries that 
perform well at off-design conditions. Further development of 
mesh movement strategies which enable large changes in the geom- 
etry are also needed in three dimensions. 

Acknowledgments 

The authors would like to thank Jamshid Samareh for his assis- 
tance in connecting the software with his geometric parameteriza- 
tion routines, and Shahyar Pirzadeh for providing numerous meshes 
for three-dimensional testing purposes. 


10 


Appendix A: Governing Equations 

Flow Equations 

The governing equations are the three-dimensional Reynolds- 
averaged Navier-Stokes equations. These equations are given as: 


. h)dQ - jkr v n)dD. 


= 0 


( 10 ) 


where h is the outward-pointing normal of the control volume 
boundary and the vector of conserved variables, Q , and inviscid 
and viscous flux vectors P t and P v , are given by 


Q = 


P 

p u 
pv 
pw 
E 


( 11 ) 


M 2 

^zz = {V- + V-t)^j\.2xv z -{u x +v y )] 


*xy ~ V " + + 


M 

Ixz = Izx = + + 


tyz ^ zy M-f) (v z Wy) 


-M x 


<lx 


J£_ + j^.^ 2 

Re(y-\)\Pr Prjdx 


-Moo f J£_ + M-t 
Re(y-l)\Pr Prjdy 


(17) 

(18) 

(19) 

( 20 ) 
( 21 ) 
( 22 ) 


Pt= 


where 


pw 


pv 


pw 

pw 2 + p 


pvw 


pww 

pwv 

i + 

pv 2 + p 

j+ 

pwv 

pww 


pvw 


pw 2 + p 

(£ + p)u 


(E + p)v_ 


( E + p)w 


-M x 


Pv = fvi + gvj + Kk 


fv 


UX xx + VT xv + WX xz - <tx 


UX xy + VX yy +WX z . 


K= 

\ux xz + vx yz + wx zz -q zA 

The shear stress and heat conduction terms are given by 


( 12 ) 


(13) 


(14a) 


(14b) 


(14c) 


~[2u x -(v y + w z )] 

(15) 

~[2v y -(u x + w z )] 

(16) 


(lz 


iL + iiiW 2 


(23) 


Re(j - 1 ) \Pr Pr t )dz 
The equations are closed with the equation of state for a perfect gas 


/ 1 \Tz 7 (w 2 + V 2 + W 2 )l 

p = Cy-l)|£-p^ -J 


(24) 


and the laminar viscosity is determined through Sutherland’s Law: 

\ 3/2 


A. = (1 +C*){T/T„y 
Aoc r/:L + c* 


(25) 


where C* = 198.6/460.0 is Sutherland’s constant divided by a 
free-stream reference temperature, which is assumed to be 460 °R . 

Turbulence Model 

For the current study, the turbulence model of Spalart and 
Allmaras 27 is used. This is a one-equation turbulence model given 
as 


n ri 

~Dt = [(v + (i + c 62 )v) v w]-c 62 \)V 2 d} (26) 

M x ( c f>, „ Yt*V , ,, , . ~~ , Re , , TJ0 

Re{ Cw ' fw k^AJLJ Cb P ft * )S 


where 


and 


/ = * 3 
U 1 Z 3 + C Vj 


t> 

X = v 


S = S + -2- — f 
S S ReK 2 d 2Jv 2 


f =1 4= 

/V2 1+Z/v, 


(27) 

(28) 

(29) 

(30) 
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In these equations, S is the magnitude of the vorticity, and d is the 
distance to the nearest wall. The function f w is given as 


fw 



1/6 


(31) 


where 

g = r + c„ 2 (r 6 - r) (32) 

and 


/ — „ i jji 

Re SK 2 d 2 

The last term in Eq. 26 is used when specifying the transition loca- 
tion. Because the computations in the present work are all assumed 
to be fully turbulent, this term is not used. Therefore, definitions of 
f t and f 2 , which are associated with these terms, are not given. 
Alter Eq. 26 is solved for t) , the eddy viscosity is computed as 

[l t = px> t = p X>f Vi (34) 

Appendix B: Linearization of Incompressible Solvers 

For both two and three dimensions, the incompressible versions 
of the flow solvers have also been differentiated. The sensitivity de- 
rivatives are verified by comparing with derivatives obtained using 
finite differences. 

For the two-dimensional incompressible flow solver, the test 
case is the 2-element airfoil shown in Fig. 1 . The angle of attack is 
1 ° and the Reynolds number is 5 million. The design variables are 
identical to those used for verifying the sensitivity derivatives for 
the compressible case. Comparisons between computed values and 
derivatives obtained using finite differences are shown in Tables 
B1-B4. Derivatives of the lift coefficient with respect to vertical 
changes in the shape design variables are shown in Table B1 with 
the corresponding derivatives for drag shown in Table B2. The de- 
rivatives obtained for translation of the flap are shown in Tables B3 
and B4. As seen in the tables, the derivatives obtained using the ad- 
joint formulation exhibit excellent agreement with the finite differ- 
ences. 

Table B1 . Accuracy of two-dimensional incompressible derivatives 
for lift coefficient (see Fig.l) 


Table B2. Accuracy of tw o-dimensional incompressible derivatives 
for drag coefficient (see Fig. 1) 



Finite difference 

Adjoint 

% Diff. 

a 

0.00807 

0.00806 

-0.124 

Point A 

0.22574 

0.22574 

0.000 

Point B 

-0.03480 

-0.03480 

0.000 

Point C 

0.16820 

0.16820 

0.000 

Point D 

0.06251 

0.06251 

0.000 


Table B3. Accuracy of two-dimensional incompressible derivatives 
of lift coefficient for flap translation 



Finite difference 

Adjoint 

% Diff. 

x-translation 

1.3358 

1.3361 

0.022 

y-translation 

-6.5537 

-6.5537 

0.000 


Table B4. Accuracy of two-dimensional incompressible derivatives 
of drag coefficient for flap translation 



Finite difference 

Adjoint 

% Diff. 

x-translation 

0.02442 

0.02447 

0.205 

y-translation 

-0.22438 

-0.22437 

-0.004 


Table B5. Three-dimensional incompressible derivatives 



Finite difference 

Adjoint 

% Diff. 

a 

-0.00307 

-0.00307 

0.000 

Twist #3 

0.00922 

0.00922 

0.000 

Shear #3 

-0.02648 

-0.02648 

0.000 

Thickness #3 

-0.04399 

-0.04400 

0.023 

Camber #4 

-1.2594 

-1.2594 

0.000 



The mesh used for verifying the three-dimensional incompress- 
ible linearizations is shown in Fig. 6. The flow conditions are an 
angle of attack of 3° and a Reynolds number of 9 million. The 
shape design variables are the same as those used in verifying the 
sensitivity derivatives for the compressible solver. Results are pre- 
sented in Table B5, and the linearizations are shown to be highly 
accurate. 
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